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Linearizing the Heisenberg equations of motion around the ground state of an interacting quantum 
many-body system, one gets a time-evolution generator in the positive cone of a real symplectic Lie 
algebra. The presence of disorder in the physical system determines a probability measure with 
support on this cone. The present paper analyzes a discrete family of such measures of exponential 
type, and does so in an attempt to capture, by a simple random matrix model, some generic 
statistical features of the characteristic frequencies of disordered bosonic quasi-particle systems. 
The level correlation functions of the said measures are shown to be those of a determinantal 
process, and the kernel of the process is expressed as a sum of bi-orthogonal polynomials. While 
the correlations in the bulk scaling limit are in accord with sine-kernel or GUE universality, at the 
low-frequency end of the spectrum an unusual type of scaling behavior is found. 



ryj 1 I. INTRODUCTION 



Perturbing the ground state of an interacting quantum many-body system and linearizing the Heisenberg equations 
of motion for the boson Fock operators, one faces the standard problem of small oscillations. Concrete examples are 
furnished by the vibrational modes of a solid, the spin waves in a magnet, the electromagnetic modes in an optical 
medium, and the oscillations of the superfluid density of a Bose-Einstein condensate. Common to these excitations is 
that they second-quantize as bosons or bosonic quasi-particles. 

Adding some amount of disorder to the system, one may ask: what are the statistical features of the excitation 
spectrum and, in particular, which of these features (if any) reflect the bosonic nature of the quasi-particle excitations? 
O , Is there some kind of universality akin to the Wigner-Dyson universality known from other disordered systems? If so, 
what are the universal laws, and what's the role of symmetry in determining these laws? 
| In the parallel case of fermionic quasi-particles the situation is now fairly well understood. If the system is of 
. metallic type and in the ergodic limit, the statistical behavior at high energies is in accord with the universal laws of 
CO Wigner-Dyson statistics. For low excitation energies, however, the canonical anti-commutation relations obeyed by 
the fermion operators make themselves felt: they constrain the form of the Hamiltonian matrix and thus give rise to 
several new universality classes beyond Dyson's threefold way pj. Some of these are realized by chiral Dirac fermions in 
a random gauge field Q , others by quasi-particles in disordered gapless superconductors 0-01 A complete symmetry 
classification of quadratic fermion Hamiltonians has been carried out p| , and the role of Riemannian symmetric spaces 
\ and superspaces in providing an effective description has been emphasized 0, Q ■ 

Progress has been slower for bosonic systems, and so for good reason, as these are set apart by several distinctive 
features from other random problems studied and solved in the past. For one thing, in the case of bosons it makes 
little sense to choose - as one often does for fermions - the matrix elements of the quasi-particle Hamiltonian as 
I 1 independent and identically distributed random variables. In fact, most of the boson Hamiltonians produced in such 
a manner would generate runaway dynamics rather than oscillatory motion around a stable ground state. In the case 
of bosons one therefore has to pay attention to the fact that the matrix elements depend in a complicated way on 
the ground state of the many-boson system and hence on the disorder of the microscopic parent problem 0, @ • As 
a technical consequence, a direct analog of the so-called Gaussian Ensembles, which were pivotal in initiating the 
Wigner-Dyson theory and establishing its universal statistics, is unavailable in the context of bosons. 

For another complication, low-frequency bosons are usually insensitive to weak disorder. Many of the excitations 
listed above are Goldstone bosons associated with a spontaneously broken symmetry, and for such excitations low 
frequency is tantamount to low wave number, or large wave length, which causes the scattering by disorder to be 
suppressed, as the disorder is effectively seen only on average over regions of size given by the large wave length. Thus 
the disorder averages out and becomes less effective, and hence the behavior of weakly disordered Goldstone bosons 
tends to be system-specific. (Of course this still leaves it possible for weakly disordered bosons of non-Goldstone type 
to exhibit universal statistics In order for any universality to set in, the disorder strength often has to be so 

large that standard calculational tools such as the impurity diagram technique fail to apply. 

In the present paper we are going to introduce and completely solve a simple random matrix model of disordered 
bosonic quasi-particles, which we believe to be most closely analogous to the Wigner-Dyson Gaussian Ensembles while 
retaining the crucial features of bosonic statistics and stability of the motion. In a follow-up paper we will investigate 
the question whether this simple model might be representative of a whole universality class of related problems. 
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To formulate the model, let qj ,pj (j — 1, . . . , N) be a canonical set of position and momentum operators, and 
consider their linearized Heisenberg equations of motion in the most general form: 



N N 



i=l 



where Bij = Bji , dj = Cji , and A y - are real numbers. If the system was invariant under time reversal (qj t— ► qj , 
Pj i— > —pj), the coefficients Ay would have to be zero, but we here consider the generic case without symmetries. The 
criterion for stability of the dynamics is that the stability matrix be positive: 

B A 



A C) >0 - 

Assuming ft* = ft > , the generator of the Heisenberg time evolution, 

_(A -B 
X - \C -A 1 

has eigenvalues that come as imaginary pairs Ouij where uij > (j = 1, . . . , N) are the characteristic frequencies (or 
single-boson energies) of the small-amplitude motion. In a classical setting one would introduce the generator X as 
the symplectic gradient of the Hamiltonian function linearized at a stable equilibrium point of the classical flow. 

The natural transformation group of the problem at hand is the real symplectic group in 2N dimensions, Sp 2 7v(B0> 
acting by linear canonical transformations on the operators qj ,pj and by conjugation on the generator X . We can 
now explain one of the distinctive features of the present problem: when formulating the Gaussian Ensembles of the 
Wigncr-Dyson theory one makes the postulate that the transformation group of the problem (Oat, Ujv, or USp 2 jv, 
as the case may be) also be the symmetry group of the chosen probability measure, whereas in our case no such 
simplification is possible. Indeed, Sp 2 jv(^) is non-compact, and a probability measure dfx cannot be invariant under 
a non-compact group action and at the same time have total mass J d/j, = 1. 

One is therefore looking for some construction principle other than symmetry. Our key here is the positivity of 
the real symmetric stability matrix ft : a natural way of building positive real symmetric matrices h is by adding a 
sufficient number of rank-one projectors with positive weights. Equivalently, we may put 

M 

hjj = VigVjg (i, j = 1, ■ • . ,2N) (1.1) 

a=l 

for some set of real numbers Vi a . We now consider the Vi a as the fundamental variables, and choose them to be 
independent and normal (or Gaussian) distributed random variables with zero mean and variance t . Then we use 
(II. 1|) to push forward the normal distribution for the Vi a to a probability distribution dfi(h) for h (and hence for X). 
If M > 27V, the result is 

dn{h) oc e-^^Det^)^'- 1 ) TT dhu , I = M — 2N > , (1.2) 

■ L± i<j 

with the domain for h still defined by h > . The probability distribution H.2fl is the object of study of this paper. 

We now give a summary of the contents and the results of the paper. After collecting some basic facts from 
symplectic linear algebra in Sect. ^ we reduce d/i(h) in Sect. II I II to a probability distribution on the space of 
characteristic frequencies lu\, . . . ,ojn (the positive eigenvalues of — iX), and find this to be 

N 

d(J,N,l(ui, ■ ■ -,uj n ) = c N j(t) - ujj){ujf - uj 2 j) Y[ t4e~™ fc duj k . (1.3) 

i<j k=l 

Using the method of bi-orthogonal polynomials we show in Sect. IV CI that the n-level correlation functions of this 
probability distribution are of determinant type and are completely determined in the usual way - see 1)5. 19[) - by a 
certain kernel Kn(uj,uj) given as a sum over bi-orthogonal polynomials. We compute the large- TV asymptotics of this 
kernel in the bulk of the spectrum (in Sect. IV D|) and at the 'hard' edge u = (Sect. IVE|) . using a contour integral 
representation of the bi-orthogonal polynomials (Sect. IVT^) . In the former case we establish the scaling limit 

sin (irp no (x)(LJ — Cj)) , w ~s 

T lim K N (Nx/T + LU , Nx/t + Q) = V , J \. ^ e -r(»)(a»-a») ) ( L4 ) 

N—too 7r(w — UJ) 
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which is independent of I . The function Poo(x) of the scaling variable x — lot/N is the large- N limit of the level 
density. Viewing ■np^x) as the imaginary part of a Green's function lim e ^o+ d( x + the function r{x) is the 
real part. We compute Poc{x) by two independent methods (from a variational calculation in Sect. IIVI and from 
bi-orthogonal polynomials in Sect. IV D|) . with the result being 

Poo{x) = ^(x/b)- 1 / 3 ^ + y/l- x 2 /6 2 ) 1/3 - (1 - - x^/b 2 ) 1 / 3 ) (0<x<b = 3V3) . (1.5) 

Apart from the last factor, which is irrelevant since it cancels on passing to the level correlation functions, the 
right-hand side of (|1.4|) is the famous sine kernel known from systems with unitary symmetry. Thus we recover 
Wigncr-Dyson universality of the class of the Gaussian Unitary Ensemble (GUE) at bulk frequencies. 
At low frequencies to ~ N~ 1 / 2 we find convergence to an unusual kind of scaling limit: 

o 2 2 

lim N- l l 2 K N {N- x ' 2 yi /T,N- l ^y 2 /T) = / du & — e -^ u+V2 / v (v/u) 1 — 9 ^ , (1.6) 

AT— too 2ir z J J v u z — v z 

iR+e Ui 

where Ui denotes the unit circle in C , and iR + e is any axis in the right half plane parallel to the imaginary axis. 
The result (ll.fcjjl is reminiscent of formulas obtained by Efctov's supcrsymmctry method, with u and v playing the 
role of radial polar coordinates of a Ricmannian symmetric supcrspace. We intend to elucidate this connection in a 
future publication. 



II. THE HAMILTONIANS OF STABLE MOTIONS 



Let there be some position variables q±, . . . , qjy and canonical momenta Pi, . . ■ , jjjy , and consider Hamiltonians H 
of the quadratic form 











-A' ) 





H = - ^2 ( C v P*Pj + B v 9iQj + A ij {liPj + PjVi)) i 
where A, B, and C are real matrices satisfying B = B l and C — C t . Rewriting H as 

H =l(l P) 

we see that the matrix, X , of H satisfies the linear condition 

x-j + jx = o, j-(j> H "J"), x-{*Z%)- M 

This is saying that X lies in sp 2 Ar(K), the Lie algebra of the real symplcctic group defined by 

Sp 2JV (M) = {g G GL 2N (R) \ g l Jg = J} . 

A matrix X £ sp 2 Ar(R) need not be diagonalizablc (e.g., the generator of free motion, A = B = and C = ljy , isn't); 
and even if it is, the eigenvalues will in general be complex. 
We now impose the condition 



B A 

A 1 C 



> , (2.3) 



i.e., we require all eigenvalues of the real symmetric matrix h to be positive. The corresponding domain in sp 2 w(K) 
will be denoted by <£°: 

€° := {X G sp 2W (M) | X = hJ , h = h l > 0} . (2.4) 

Although the eigenvalues of h have no direct relation to the dynamics of the system, positivity of h ensures that the 
motion generated by the Hamiltonian H is stable, or 'elliptic'. As a consequence of ellipticity, there exists some linear 
canonical transformation (q , p) — > (Q , P) which takes the Hamiltonian to a sum of harmonic oscillators, 

1 N 

i=l 
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with ojf > . Put differently, for X £ €° one can always find a symplectic transformation g £ Sp 2JV (R) that conjugates 
X to quasi-diagonal form: 

X = gflg' 1 , °=f^ 0^) ' u = diag(wi,w 2 , ■ ■ ■ ,wn) , (2.5) 

with real and positive cjj (i = 1, . . . , iV). 

All of the discussion below will be based on the elliptic domain (E . Let us therefore collect some of its mathematical 
properties. First of all, if X is in £°, then so is its conjugate gXg~ Y by any element g £ Sp 2JV (IR). Thus €° is invariant 
under the action of Sp 2Af (IR) on <E by conjugation. Secondly, let t denote the Abelian algebra of block-diagonal 
matrices of the form of £1 in i|2.5[) but with diagonal elements Wi that are any real numbers (not necessarily positive). 
Let t+ C t be the subset of block-diagonal Vl with positive u>i . Then, as we said earlier, every X £ €° is conjugate to 
a unique Q £ t+ by some g £ Sp 2Ar (IR). Thirdly, introducing T := exp(t), which is an AT-dimensional compact torus, 
T = (S 1 )^, let G/T be the quotient of G = Sp 2Ar (R) by the right action of T. Then the mapping 

(G/T) x t+ £° , (gT, ft) ^ gtlg- 1 (2.6) 

(the reverse of the process of quasi-diagonalization), is a smooth bijection. 
We are stating these facts without proof, as they are standard facts of symplectic linear algebra. 



III. PROBABILITY MEASURE 

By placing a probability distribution on the elliptic domain (£°, one gets a random matrix model for disordered 
bosonic quasi-particles. We are then interested in the statistics of the characteristic frequencies or levels u>i . 

It is well known that in the Wigner-Dyson situation of random Hcrmitian or random real symmetric matrices, 
where the symmetry group is compact, the level correlation functions exhibit universal behavior in a suitable scaling 
limit. One may therefore ask whether a similar scenario - leading to universal laws, possibly of a new kind - might 
be at work in the case being considered. 

To answer this question we need to investigate a class of probability distributions on (5 as wide as possible. As a 
first step, the present paper deals with a family of well motivated distributions which are easy to analyze. 



A. Choice of measure 



Coming from the standard Wigner-Dyson situation with a compact symmetry group, one might be inclined to try 
and consider a Gaussian distribution 

P(X)dX oc c'^^dX , 

where dX is a Lebesgue measure for £°: 

dX := Y[ dAij . dBij dC l3 . (3.1) 

However, such a distribution has infinite mass, since it is invariant under the action X i— > gXg~~ x by the non-compact 
group Sp 2JV (]R), and it therefore cannot be normalized to be a probability measure. 

Staying within the class of Gaussian distributions, a better choice of distribution function is 

P(X = Jh) CX e-^Va-'T"' = c -rTr(J^X)/2-,Tr(J^Xf (3 2) 

for some positive parameters cr,r. Because of the presence of J -1 under the trace, this distribution function is 
invariant under conjugation X t—. ► gXg~ x only if g £ Sp 2W (R) satisfies the additional condition g~ x J g = J. Combining 
the two conditions, g t Jg = J = g~ 1 Jg , one sees that the invariance group of the function P(X) in l|3.2|l is the 
intersection of the real symplectic and orthogonal groups in 27V dimensions: 

K = Sp 2W (M) n S0 2A r(K) . (3.3) 

This group K is isomorphic to Uat, the group of unitary transformations in N complex dimensions. Indeed, changing 
from the symplectic basis {qi, . . . , qN,Pi, ■ ■ ■ ,Pn} to the oscillator basis 

{ai,...,a N ,a\,...,a f N } , a,j = -j={qj + ipj) , a] = -^=(qj - ipj) , 
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one finds that K is the subgroup of canonical transformations that do not mix the lowering operators {a^} with the 
raising operators {a]}- Moreover, XJn — K C Sp 2Ar (]R) is known to be a maximal compact subgroup. It therefore is 
the biggest symmetry group possible in our problem. 

In the sequel we will consider (|3.2|) with a = . Thus we take our probability distribution to be 

P(X)dX :=c N (T)e- rTl{J ~ lx)/2 dX , (3.4) 

with the normalization constant cjv(t) chosen in such a way that Jg P(X) dX = 1. Further motivation for this choice 
of distribution was put forth in the introduction (Sect. 01. 



B. Polar decomposition and reduction 

Let now F(X) = F(gXg^ 1 ) be some function on €° which is radial, i.e., invariant under conjugation by every 
element g € Sp 2 n(^)- Given such a function F, which depends only on the eigenfrequencies uj\, . . . , ojn of X, we wish 
to take the expectation of F w.r.t. the probability measure P(X) dX : 

(F) := [ F(X)P(X)dX . (3.5) 

The problem of computing such expectations is best tackled by using the polar decomposition (£° = t+ x (G/T) which 
is given by quasi-diagonalization of X; see (|2.6|) . Inserting that decomposition into l|3.5|l one has 



(F) = J PigVLg' 1 ) dg^j F(Q.) j(Q) dfi , rfO = dwxdw 2 ■ ■ ■ duj N , (3.6) 

where gr is a G-invariant measure for G/T, and is the Jacobian of the change of variables X = gtlg -1 . 
Let us calculate this Jacobian. Differentiating the polar coordinate mapping l|2.6() we get 

6( g ng- 1 )=g(6n+[g- 1 6g,n})g- 1 ■ 

The Jacobian we are seeking is the product of all nonzero eigenvalues of the linear operator X i— > [X, £1} . These 
eigenvalues are called the roots of the pair (sp 2 jv W, t). They are 

±(u) t + UJ 3 ) (i < j) , ±((J i - ujj) (i < j) , 

each with multiplicity one. Thus, by taking the product of all non-vanishing roots, 

TV 

i(O) dil = \{{ujj - uj]) 2 \{(2uj k ) 2 dw k . (3.7) 

i<j k—1 

To complete the polar integration formula l|3.6|l we need J G ^ T P(gflg~ 1 ) dgx- In the next subsection we are going 
to show that this integral can be calculated in closed form and depends onwi . . . , ojn as 

r N 

/ Pigttg- 1 ) dg T oc T\(ui + (Jj)- 1 TT 1 e"™ fc . (3.8) 
Thus, in total, the expectation of a radial observable F(X) = F(Q,) = F{u>\, . . . ,u>n) becomes 

N 

(F) = cjv(t) / n F(u u ...,w N ) Ufa - uj,)^ - w?) e"™^ dw k , (3.9) 

•^ R + i<j k=l 

where cat(t) is another normalization constant. This expectation, for the special choices of F that give the level 
correlation functions, will be calculated in Sect. El of the paper. 
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C. Computation of the integral 113.811 

We now establish (|3.8() . Omitting a normalization constant, we denote the integral on the left-hand side of l|3.8[) by 

7(0):=/ e-^^v^^dgr . (3.10) 
Jg/t 

What makes this integral computable in closed form is that J lies in sp 2 jv(^) an d O i— > pfigr -1 is the adjoint action of 
G = Sp 2 7v(^) on its Lie algebra. These circumstances place the integral in the class of integrals of Harish-Chandra- 
Itzykson-Zuber type, which are covered by the Duistermaat-Heckman theorem and its generalizations. In the present 
case, the integral can be computed in a particularly simple manner, as follows. 

Let dg and dt be Haar measures for G and T, respectively, with dg = dgT dt and J T dt = vol(T). Our first step is 
to switch from G/T to integrating over the full symplectic group G : 



7(0) = _L_ f e -TTr(j-i g n g -i)/2 dg 
'g 



vol(T) 



Next we use that dg is invariant under inversion, g > g^ 1 . After this transformation the integrand is expressed in 
terms of the combination gj~ 1 g~ 1 = —gjg~ 1 . Since kJk^ 1 = J for k e K = \Jn , we can push down the resulting 
integral over G to an integral over the quotient space G/K. Let dgx and dk be invariant resp. Haar measures for 
G/K and K so that dg = dgx dk . Then 

/(0) = ^ffl/ e^W'-^dgK, vol(K)=[dk. (3.11) 

V0l(T) J g/k Jk 

The homogeneous space G/K = Sp 2 jv(M)/Uiv has the salient feature of being a non-compact symmetric space 
of Hcrmitian type. Such spaces carry the structure of a Kahler manifold, which means that G/K comes with a 
non-degenerate, closed, and G-invariant 2-form (the Kahler form of G/K). Writing gJg~ x =: Q this is the form 

f3 = Tr (QdQ A dQ) . (3.12) 

Notice that dim R G/K = N(2N + 1) - N 2 = N(N + 1). Raising f3 to its ±N(N + l) th exterior power one obtains a 
top-dimensional form, pi N ( N+1 \ which is still G-invariant and nonzero. Since G/K is homogeneous, there can be at 
most one such form up to multiplication by scalars. Therefore, there exists some (nonzero) constant such that 

dg K = const x / 33 JV ( W + 1 ) . (3.13) 

By Darboux's theorem one can find local symplectic coordinates for G/K that bring j3 into canonical form. While 
this fact by itself would not be of much practical help, in the present case such coordinates exist globally and, moreover, 
they can be chosen in such a way that Tr(0 gjg -1 ) depends on them quadratically. 

To describe these perfect coordinates, consider the space of complex symmetric N x N matrices, Sym(C Ar ), which 
has dimension iiV(iV-l-l) over C and thus shares with G/K the dimension N(N+1) over K . With every Z € Sym(C Af ) 
associate a positive Hermitian 2N x 2N matrix g by 

^ S(Z , Z ^(^ Z z f' 2 (1 + z V). 

Now if S is the matrix of the unitary transformation from the real symplectic basis {pj , qj } of IR 2 ^ to the oscillator 
basis {oj , a]}: 

g 1 ( i In 



^2 V -1 ^ i!jv, 

then g := S~ 1 gS is immediately seen to be a real matrix and, using the relation 



SJS- 1 = iS, , E, = 



In 
-1 N 
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one finds that g = S 1 gS satisfies Jg = g l Jg = J and hence lies in Sp 2 jv(K). Moreover, the reverse correspondence 
k h- > k = SkS -1 is the isomorphism between K and U_/v discussed in the paragraph after 1)3.3(1 : it takes k G K to the 
block-diagonal form 

* = ( o v) ' u e Vn • 

It is now clear that the mapping Sym(C Ar ) — > G/7£T by Z i— * S^ 1 g(Z, Z^)S ■ K = gK is a bijection. Using it to express 
the Kahler form (3 in terms of the complex symmetric matrix one obtains 

(3 = Tr (Q dg A dQ) = -4i Tr (E 3 d<T x A dg) = 8i Tr (dZ A di^ ) . (3.15) 

Thus, the top-dimensional form fih N ( N + 1 ) is constant in Z : 



_ (8i) ^(iv+i) 2 4JV(Ar-i) tt dz d2 (3 16) 

(|JV(JV + 1))! iJ -'<j J J 

and from (|3.13|) the invariant measure dgx is a constant multiple of the Lebesgue measure for Sym(C Ar ). 
Finally, from g = S^^^gS , SJS^ 1 = iS 3 and l|3.14|l one has 

-Tr(fi ff J 5 " 1 ) = Trw(l + 2ZZ t ) + Trw(l + 2Z t Z) , u = diag(wi, . . . , w N ) . 

Our integral 1(3.11(1 now becomes a Gaussian integral: 

7(0) = const xe"^' u ' / e -37>* rftfu+uZ) JJ dz .. ^ , 
Doing this integral one immediately obtains the result for 7(0) stated in 13.8(1 . 



D. Generalization 

A slight generalization of 1(3.4(1 is afforded by the observation that the determinant of X in ((2.2(1 is always positive: 

N 

Det(X) = Det(O) = u\ . 

fe=i 

Thus, by multiplying the probability measure P{X) dX by some power I— 1 > — 1 of the positive square root Det(X) 1 / 2 
and adjusting the normalization constant, we get another probability measure: 

Pi{X)dX = const x Dct(X)^ 1 - 1 ^ e-^ TTl(J ~ lx) dX . (3.17) 

This measure is still Uiv-invariant. By the process of quasi-diagonalization and drawing on our results above, we push 
it forward to a measure for the eigenfrcquencies. The result is 

N 

dn Nj i(u>\, . . . ,uj n ) = c Nj i(t) - ujj) 2 ^ + uij) Y[ e _ ™ fe duj k . (3.18) 

Kj fc=l 

This, for any non-negative power I <G Z , is the family of probability distributions to be studied in the present paper. 



IV. LARGE- LIMIT OF THE 1-POINT DENSITY IN THE BULK 

The 1-point density p(lu) dio is defined as the probability density for any one of the eigenfrcquencies to have the 
value of uj . irrespective of what the values of the other eigenfrcquencies are; thus p(iv) is the function 

N 

p(u>) ■■= / y]S(u; - Wj) dfi N) i(u;i, . . . , (4.1) 
^ i=l 
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FIG. 1: The graph of the function p a 



which has the properties p{u>) > and 



p(to) dcu = N 



(4.2) 



We are now interested in the behavior of the density function p{uj) in the limit of N — > oo. From the expression 
(|3.18jl and experience with similar problems (see, e.g., UK), we expect that this limit can be obtained by maximizing 
the functional 



1 

F = - 



OO pOO 



\n{{uj - lo') 2 {uj + lo')) p(uj)p(uj')duj'duj + / ln(u/ e"" T ) p(uj) du 



(4.3) 



subject to the constraint (|4.2(l and the condition p(uj) > 0. More precisely, the limit is expected to exist in the scaled 
variable x := ujt/N; i.e., there should exist a certain non-negative function Poo{x) with J Poo(x) dx = 1 such that 
p(uj) is asymptotic to t poa(ujT / N) . 

Varying F with respect to p(to) we get 



5F 



(2 In \uj — cj'| + ln(cj + cj')J p(w') rfw' + Z In w — 



We now insert the asymptotic equality p(Nx/r) sa Tp oc (x) and pass to the limit iV — > 00 in the scaling variable x . 
Let supp(poo) = [0,6] be the region of support of poo ■ Then the condition 5F/5p(uj) — NX, where A is a Lagrange 
multiplier for the constraint l|4.2|) . yields the equation 



(2 In \x — x'\ + ln(x + x')) Poc(x') dx' — x = X (0 < x < b) , 



(4.4) 



which no longer depends on the parameter / . It can be shown that our functional F is convex; as a result, the solution 
Poo of Eq. H4.4f) exists and is unique when supplemented by the normalization condition 



Poo(x) dx = 1 



(4.5) 



In the following subsections, we are going to construct the solution to the mathematical problem posed by l|4.4() 
and l|4.5|) . It will turn out to be 

Poc(x) - -^(a;/6)- 1/3 ((l + y/l- x 2 /6 2 ) 1/3 -(l-^l- x 2 /6 2 ) 1/3 ) (0 < x < b = 3^3) . (4.6) 
The graph of this function is plotted in Fig. ^ From the expression (|4.6() the behavior near the lower edge x = is 



Poo (x) ~ ^(26/z) 1/3 



(0 < x < 6) 
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while close to the upper edge x = b one gets 

Poo (x)~^-(l~x 2 /b 2 ) 1/2 (x<b,x^b). 

In the vicinity of the upper and lower edges there exists crossover to a fine-scale behavior that cannot be found by 
the present method of maximization of the functional F. The crossover at the upper edge involves Airy functions on 
a scale iV 1 / 3 , which is small compared to the bulk scale N. At the lower edge, the crossover occurs on a very fine 
scale, iV -1 / 2 , which is small even in comparison with the bulk mean level spacing (which is of order N°). 



A. Method of solution (idea) 

We do not know how to solve Eq. (|4.4|) for the unknown function Poo{x) directly. Therefore, to simplify the problem 
we differentiate once with respect to x to obtain the equation 

2V f b Poojx^dx 1 + r b Poc (x')dx' =1 

Jq X — X 1 Jq x + x 1 

where V means the principal value of the integral. At this stage, the value of b is unknown but assumed to be finite. 
Introducing the Green's function (or Stieltjes transform) 

9(z):= fpM^ 7 zeC\[0,6], (4.8) 
Jo x - z 

and the related functions 

g±(x) := lim g(x±ie) , g Q (x) := -g(-x) , (4.9) 

we bring Eq. (|4.7|l into the form 

g+(x) + g.[x) + g {x) = -l (0 < x < b) . (4.10) 

To solve this equation, we are led to do an exercise in complex analysis which is motivated as follows. 

Let w i— > f(w) be some meromorphic function of a complex variable w, and let the equation z = f(w) have r simple 
roots w\{z) ,W2{z) , . . . , w r (z), i.e., z = f(w\(z)) = . .. = f{w r {z)). If the function / is analytic in 1/w at w = oo, 
these roots add up to a constant: 

r 

Wj (z) = const =: c (independent of z) . (4-11) 

i=l 

Indeed, if 7 is a closed contour encircling all of the roots in the counterclockwise sense, then 

En \ - 1 1 / dw 

where the second equality is by the residue theorem, and the last equality follows by contracting 7 to the point at 
infinity. Thus ^w'^z) — and hence ^Wi(z) = const. 

Eq. (|4.11|) for r — 3 looks similar to (|4.1U|) and can, in fact, be made to look identical to it by the following 
observation. Notice that the function z 1— > g(z) defined by (|4.8fl is holomorphic in the interior of the left half of the 
complex plane. Suppose, therefore, that we have found a root g{z) of z = f(g{z)) which is holomorphic in the left 
half plane, and that g±{x) = lim e _»o+ g(x±ie) are its two analytic continuations to positive real x £ (0, b). Moreover, 
suppose that the function / has a reflection symmetry 

f(w) = -f(2a-w) (oeC). (4.12) 

Then z 1— * 2a — g(—z) is a root of z = f(w) holomorphic in the right half plane, and from l|4.11l) we infer that 

9+(x)+g-(x) + (2a-g(-x)) = c . 

Setting go(x) = —g(—x) this becomes the same as (|4.10|) if 

c-2a=-\. (4.13) 

Thus we are inspired to interpret g + , g_ , and 2a + go as the three roots of an equation z = f(w). Given this 
interpretation, solving l|4.10|l amounts to finding the function /. 
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B. The good function / to consider 



We are looking for a certain meromorphic function / on C . By adding a point at infinity we can view such a 
function / as a mapping of the Ricmann sphere S 2 = C U {00} to itself. We want this mapping to have degree r = 3 ; 
i.e., every regular point z of / is to have three distinct pre-images: 

f~ l { z ) = {■wi(z),w 2 {z),w 3 (z)} . 

Such a mapping can be presented in the general form 

f(w) = foo + J2 ( 4 - 14 ) 

with some complex numbers a; , 6; , and /oo . 

Let us narrow down the choice of parameters. From the normalization condition (|4.5J) and the definition of g(z) in 
(|4.8|) . we have the limit zg(z) — + — 1 for z — > 00 . Therefore, since f(g(z)) = z by construction, we need f(w) to have 
a pole at w = g(oo) = with residue —1. So we choose a\ = and b\ = —1. The reflection symmetry (|4.12|> is then 
implemented by setting /oo = , 63 = b\ , and a% = (i — 1) a for i = 1, 2, 3 and some a € C . Thus, 

/ H = + 7T > 

u> w — a w — la 

where the parameters a and 62 are still unknown. 

Next observe that for a degree-r holomorphic mapping / : S 2 — > S 2 , the number of singular points, where f'(w) = , 
is 2r — 2 . Indeed, writing / as f(w) = p{w)/q(w) where p and q are polynomials of degree r , one has 

ti, x p'(w)q(w) -p(w)q'(w) 
J [w) = ^ - 

the numerator of which is a polynomial of degree 2r — 2 and so has 2r — 2 zeroes. 

Thus we should expect our function (|4.14(l to have 2x3 — 2 = 4 singular points. The reflection symmetry (14.121) 
makes for their images {f(w) £ C| /'(?«) = 0} to be arranged symmetrically around z — 0. Now notice that our 
Green's function g(z), being the Stieltjes transform of Poo(x) with support [0, b] , must have singularities at z — and 
z = b. The image of the singular set had better contain these values, and thus is determined to be {—6,0, +6} by 
reflection symmetry. Actually, since our situation calls for / to have four singular points, the singularity at z = 
(corresponding to w = 00) must have multiplicity two. This is achieved by choosing b 2 = —61 — 63 = +2 , so that 

12 1 -2« 2 
f(w) = + — = — — , 

w w — a w — la w(w — aj{w — la) 

resulting in the behavior f(w) ~ w~ 3 for w — > 00. The singular points of / now arc w — a ± a/\/i , and 00. These 
correspond to z = f(w) = ±3-\/3/a , and , respectively, so we infer 

6 = 3V3/|o|. (4.15) 

It remains to pin down the last unknown parameter a . For that purpose, recall that the sum of the roots f~ 1 {z) — 
{wi(z), W2(z), W3(z)} is a constant, c, independent of z. To determine this constant, look at J2 w i{°°) and use that 
the poles of / are at w = , a , 2a to obtain 

c = Wj{z) = Wj(oo) = 3a . (4.16) 

We then conclude a = — 1 from P~T3)l . and b = 3^3 from jPfy In summary, the good meromorphic function / for 
us to consider is 

w " f ^ = ( ZT^ ^ ' ( 4 ' 1? ) 
w(w + l)[w + 2) 

Let us mention in passing that the idea to consider the equation z = f(w) or, equivalcntly, 

w(w + l)(w + 2) + 2/z = , 

first came to one of us (HJS) from previous work [Icj on the Green's function of the Bures measure, whose large-iV 
limit leads to a similar equation. 
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C. Solution of the problem 

The situation can now be succinctly described like this: thinking of 

W :=C\ {-1 + 1/ Vs, -1-1/^3} , Z:= (C \ {b,0,-b}) U {00} (6 = 3^3), 

as two Riemann surfaces W and Z, the function / of l|4.17|l gives us a holomorphic cover 

f:W^Z, f-\z) = { Wl (z), w 2 (z), w 3 (z)} . 

What's the monodromy of this cover, i.e., what happens when the locally defined functions z 1— > w%{z) are analytically 
continued around one of the singular points z — b ,0 , ~ b? At the point z — (or w — 00) we have a cubic singularity 
z ~ w~ 3 . Consequently, the monodromy at z — cyclically permutes the roots Wi(z). Turning to z = ±0, we see 
that linearization z = ±0 + Sz and w = f (±6) + Sw gives 

<fe - (8w) 2 . 

In the latter two cases the monodromy must exchange two of the Wi{z) while leaving the third one invariant. 

Now focus on the situation near the singular point z = — and denote by w{z) = g(z) the root which, there, 
is trivial under monodromy and hence exists as a holomorphic function in some neighborhood of z = —b . With 
the remaining two singularities being at z = and z = 6, the function g(z) actually extends to a holomorphic 
function on the Riemann sphere C U {00} cut along, say, [0, 6] CR. Let us verify that this holomorphic function 
g : (C \ [0, b}) U {00} — > W coincides with the Green's function l|4.8[) solving our problem (14.10(1 . 

By the holomorphic nature of g and Cauchy's theorem, we have that 

If g{z')dz' 
27ri J 1 z' — z 

where 7 is a small loop running around z in the counterclockwise sense. Since g is holomorphic at infinity, the loop 7 
can be deformed (through infinity) to a loop encircling the cut [0, b] , but now with the orientation reversed. Collapse 
the deformed loop to the two line segments connecting with b. Then, setting g±(x) = lim e _^o+ g(x i ie) and 

Poo (a) := 9+{x) ~ 9 - (x) (0 <*<&), (4.18) 

Z7T1 

g(z) is obviously given by the integral in l]4.8[l. 

Because g+(x) and g-(x) arise by analytic continuation from g(z) G / _1 (z), these are two of the three elements in 
the set f~ 1 (x). How is the third element of f~ 1 (x) related to g{z)l To see that, recall a = — 1 and from l|4.12ll the 
invariance of the equation z = f{w) under (z, if) 1— > (—2, 2a — w). Thus, if g{— z) is a root over — z , then —2 — g(—z) 
is a root over z , and it follows that go(x) — 2 with go(x) := — g(— x) (for < x < b) is a root over x . The roots g + (x), 
g-(x), and go(x) — 2 all are different as functions. In fact, 3mg + (x) > = Umpo( a; ) > 3 m 9-i x ) for < x < 6. So, 

Z -1 ^) = {g+i x ) 1 9-( x ) > 9o(x) - 2} , 

and from (|4.16|) we deduce that 

g+{x) + g-(x) + go(x) — 2 = w = c = 3a = —3 (0 < x < b) , 

which agrees with (|4.1U|) . Recall that in order for our analysis to work out we must choose 

6 = 3V3. (4.19) 

With a full understanding of the situation in hand, it is now an easy exercise to obtain poo(x) in explicit form. 
Solving the equation z — f(w) one finds the holomorphic function g{z) in the interval —b < x < to be 

g(x) = {-x)~^ 3 {l + y/l -x 2 /b 2 ) 1/3 + (-xy 1/3 {l - ^/l-x 2 /b 2 ) 1/3 - 1 , 

where all square roots and cubic roots are understood to be positive. This function indeed extends holomorphically 
to a neighborhood of x — — b , as the Taylor expansion at x — — b contains only even powers of y/l — x 2 /b' 2 . Analytic 
continuation around the singularity at z — gives 

g±(x) = x- 1/3 e ±7ri / 3 (l + y/l-x 2 /^) 1 / 3 + ar 1 / 3 e =F,ri/8 (l - y/l - x 2 /b 2 ) 1 / 3 - 1 (0 < x < b) . 



12 



Computing the difference (|4.18(l we then get the result for Poo(x) claimed in 1|4.6[) . with the value for b given by i|4.19[l . 
As a final remark, let us note that the good form of g(z) to use near infinity is 

1/3 / , ± \ 1/3 



From this, all moments of Poc(x) dx can be found by expanding g(z) in powers of 1/z. 

V. EXACT SOLUTION USING BI-ORTHOGONAL POLYNOMIALS 

We now express the probability measure i|3.18fl as 

N 

d(iN t t(ui, . . .,ujn) = cn,i(t) - - UJ?) Y[ c~ TUJk uj l k duj k , (5.1) 

i<j k—l 

and embark on another approach to handling it. 

To get started, recall the formula for the Vandermondc determinant: 



~\[u)i-u}j) = Det(w* 

i>j 

Using it we reorganize the probability measure 1)5. lfl as 



1 1 ... 1 

LJ\ Ul2 ■ ■ ■ 



N-l N-l N-l 
LOi Wn . . . UJ 



N 



dfx N ,l(^ ...,u N ) = c nj {t) Det (w*- 1 ) Det (wf" 2 ) e"™* J k dio k . (5.2) 

fe=i 

We also simplify our notation by setting r = 1. 

By standard properties of the determinant, Det (w* _1 ) changes only by a multiplicative constant when the monomials 
Ljj~ are replaced by any polynomials in uij of degree i— 1 . We have two Vandcrmonde determinants, O^jl^i ~~ 10 j) 
and n^jl^f — so wc introduce two sets of polynomials, denoting those of the first set by Pi^i(uij) and those of 
the second one by Qi-i(tu 2 ). Our measure then becomes 

N 

dfi N j(u;i, ...,w N ) = c N jDet (P 1 ^ 1 {lj j )) Det (Q^oj 2 )) J[ c~" k J k du k . (5.3) 

fc=i 

In order for the introduction of the polynomials P n (uj) and Q n (uj 2 ) to be useful we require them to be orthogonal 
with respect to the integration measure e~ w a/ dio : 

poo 

I m ,n= / P m {^)Qn{^ 2 )e^ui l duj = h n 5 m , n , (5.4) 
Jo 

where the numbers h n = depend on the choice of normalization for P n (ui) and Q n {uj 2 ). Such polynomials are 
constructed by a variant of the Gram-Schmidt algorithm, as follows. 

A. Bi-orthogonal polynomials 

We review the construction in the general setting of two real vector spaces V, W with a pairing (or non-degenerate 
bilinear form) 

7 : V x W —> R . 
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Given some basis Vo, l>i, V2, ■ ■ ■ of V, and a basis wo,w\, u>2, ■ ■ ■ of W, let the entries of the corresponding pairing 
matrix be denoted by 

7ro,n := l(v m , w n ) (m, n = 0, 1, 2, . . .) . 
The goal now is to construct a new basis eo, e\, 62, . . . of V, and a new basis /o, /1, /2, ■ ■ • of W such that 

77 — 1 77— 1 

^77 ~l~ ^ ^ V n ' , t /Vt l/J n 7 ^ ^ B nn ' W n ' , 

n'=0 n'=0 

(with real coefficients A nn i and -B nra ' ) , and the transformed basis vectors form a bi-orthogonal system: 

7( e m,/n)=0 m^n. 

This problem has a unique solution by the process of Gram-Schmidt orthogonalization. A nice way of presenting 
the solution is by means of the following determinants (where, by a slight abuse of notation, the matrix entries in the 
last column resp. last row are vectors, whereas all of the other matrix entries are numbers): 



D 



n-l 



70,0 
71,0 



7n,0 

with normalization factor 



70,77-1 V() 
7l,n-l VI 



7^,71—1 ^71 



In = D-\ 



D n = 



7o,o 



7o,o 7o,i 
771-1.0 771-1.1 

Wq Wi 



To, 



70,77 



7ti-i, 



(n=l,2,...), (5.5) 



7",0 ■ ■ ■ 777,77 

These formulas are easily verified. Indeed, pairing e m with w„ for m > n one gets 



7(e m ,w n ) = 



7o,o 
7i,o 

7777,0 



70,777-1 70,71 

71,777-1 7l, 71 



-1 7ti 



0, 



which vanishes because the last column coincides with one of the other columns. Since /„ is a linear combination of 
the vectors w n i with n' < n , it follows that 7(e m , /„) = for m > n. The same conclusion for m < n follows by 
reversing the roles of V and W. Notice that e„ = 1 • v n + . . . and /„ = 1 • w n + . . . by insertion of the factor D~^_ x . 
The non-vanishing pairing matrix elements for n > 1 are 7„. n = 7(e„ , /„) = 7(e„ , w n ) = D n / D n -i . 

To apply these general formulas to the case under consideration, we choose the vectors v rn and w n to be the functions 
co 1— ► to" 1 resp. co i— ► w 2n , and take the pairing to be given by integration with our measure e -UJ u/ dui : 



In 



m~\-2n — lj I 



e" w w' dto = r(m + 2n + Z + 1) 



(5.6) 



Making the identification e n = P n (ui), the general formula for e n in (|5.5|) then gives Po(^) = 1 and 

(n > 1) . 



T(l + 1) 

r(i + 2) 



T(l + 2n- 1) w° 
r(Z + 2n) 



r(/ + n + 1) ... T(l + 3n- 1) w 
Similarly, identifying /„ = Q ra (u; 2 ) we obtain Qo(^ 2 ) = 1 an( i 

r(z + i) r(z + 3) ... r(/ + 2n + i) 



(5.7) 



Q„(o, 2 ) = D~\ 



T(l + n) T(l + n + 2) 



T(l + 3n) 

, ,277 



(n>l) 



(5.8) 
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~[2 k k\(2k + iy. 



Using the relation T(z + 1) = zT(z) an easy Gauss elimination process gives the normalization constant as 

r(Z + l) ... T(l + 2n+l) 
D n = : •.. : 

T(l + n + l) ... T(l + 3n+l) 

From this, note the diagonal pairing matrix elements ho = l\ and 

PnM Qn{u 2 ) e^Lu 1 diu = h n = D n /D n _ 1 = 2 n n\ (2n + l)\ (n > 1) 



k=0 



(5.9) 



(5.10) 



B. n-level correlation functions 



The n-level correlation function R n (uj\, . . . ,uj n ) in the present context is defined as 

R n ((j 1 ,ui2,...,u) n )=n\ ^ $(^1 - ^ii)S{u>2 - ■■■ 8(u n - u) in )dfj, N ,i{u>i,-- -,& N ) . (5.11) 



il<i2<--<in 



A closed-form expression for it can be given from the bi-orthogonal polynomials P n >(u>) and Q n i(uj 2 ) for < n' < N. 
The result will take its most succinct form when expressed in terms of the modified functions 



P„(w) := (-2)-™7i!- 1 e^P Il ( W ) , 
Q„(w) := (-l) n (2n + l)r 1 Lu l Q n (io 2 ) , 



(5.12) 
(5.13) 



(the motivation for the sign (— 1)" will become clear later), which from (|5.4(l and (|5.10|) obey the orthogonality 
relations 



P m (u) Q n (uS) duj = 5„ 



The probability measure l|5.3[) expressed by these functions takes the form 

d/ijv,i(wi, ■ ■ ■ , u N ) = -^y Det (P i _i(w J )) Dct (Q i - 1 (uj j )) JJ fc du k 
Now, by using the multiplicative property of the determinant, we can also write 



(5.14) 



dfj, N ,i(ui, . . . ,uj n ) = ^j Dct (^( w <> w 3'))<,j=i,...,iv II fe dwfe ' 
where the kernel K{uji ,ujj) is defined by 

N-l 

K N (uii,L)j) = ^ Pn(Ui)Qn(Uj) ■ 

From the orthogonality relations l|5.14|l this kernel has the reproducing property 

/>oo 

/ K N (uji , u>) K N (u, uij) duj = K N (uji , 

JQ 

and the trace 

/*OC 

Kpf(uj, uj)duj = N . 



(5.15) 



(5.16) 



(5.17) 



(5.18) 



To proceed further, take notice of the relation 



K N (u n ,ut) ... K N (lj„ , uj n ) 



duj n = (N - n + 1) 



K N {ui,Ul) ... K N (wi, CJn-l) 
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which is proved by expanding the determinant with respect to the last row or column and exploiting the properties 



(|5.17|) and (|5.18|) . Using it, an inductive procedure starting from Rjy(u>i, . . . , u>jv) = Det (-K"jv(a>j j _ 1 N gives 



the n-level correlation functions as 



Rn(wi, ■ ■ .,w n ) = Det (K N (u) i ,u J )) t n 



(5.19) 



Thus the correlations are those of a determinantal process and are completely determined by the kernel Kjss{u>i ,u>j). 
The remaining discussion therefore focuses on this kernel, but first we make another preparatory step. 



C. Contour integral representation 



We are now going to show that the functions P n {ui) and Q n (u>) have expressions as complex contour integrals: 



S e (0) 



(1 - u- 2 )- n -W- 2 du/m 



(1 - v- 2 ) n v~ l - x dv/2m 



(5.20) 
(5.21) 



Both integrals are over circles in the complex plane with radius e and counterclockwise orientation; the first circle is 
centered at u = 1 and has radius e < 2 (to avoid the singularity at u = —1), the second one is centered at v — . 

Our proof of these expressions for P n {ui) and Q n {oj) will be indirect and in two steps. First, we establish some 
information on power series. In the case of Q n (uj) we insert the power series of the exponential function e" 11 , use the 
binomial expansion of (1 — v~ 2 ) n , and compute a residue to obtain 



?»M = £ 



fe=0 



kj (2k + l)\ 



In the case of P n (co), calculating the residue at u = 1 we have that 



2 d n ( e-^u 71 ^- 1 



\du n \{l + u- 1 ) 



t\n+l 



(5.22) 



(5.23) 



In both cases, defining P n (u>) and Q n (uj 2 ) by the reverse of the relations l|5.12l l5~T3*|) . we see from (|5.22l I5~2*3*)) that 
these are polynomials of degree n in u> resp. ui 2 and that the highest-degree term (oj n resp. u> 2n ) has coefficient one. 

Recall now from Sect. IV M that, given these properties, the polynomials P n (uj) and Q n (u> 2 ) are completely deter- 
mined by the orthogonality relations (|5.4(l for m ^ n. Via H5.12I5.13|I the latter are in one-to-one correspondence with 
the orthogonality relations l|5.14|) (still for m ^ n). Therefore, defining 



Pm 

(w) Qn(w) du , 

JO 



(5.24) 



the second and final step of our proof is to show that I m ,n = for m =/= n. 

To that end, we insert the expressions lf5~27H l5"2T|) into lpT2l)l . The w-dependence then is e~^ u - v ) with u £ S e (l) 
and v £ S e (0). Taking the radius e to be very small (e 1), we have that q-^( u - v ) decreases rapidly as u) goes to +oo. 
Therefore, the integral over ui exists, and we may interchange the order of integrations. Doing first the cj-integral, 

OO i 

u-v 
the remaining contour integrals for I m ^ n defined by (15.241) arc 

f u 1 - 2 ( f (l-v- 2 ) n dv\ du 

m ' n ~ Id, {i-u- 2 r + i [f sm v i+i( V -u) ) ^ ■ 

To simplify the inner integral over v we use the identity 

1-v- 2 ' 



l-u- 



2 o n— 1 / _o \ k 



Tl-1 

i- V(i 
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Inserting this into the expression for 7 mj „ we see that the terms in the fc-sum do not contribute as the residue at v — 
vanishes for all of those terms. Doing the w-integral for the first term on the right-hand side, we get 



so the remaining u-integral is 



v~ l ~ x {v - u^dv = -2mu~ l ~ x 



I m< n = (^r 1 f (1 - u-^-^u^du . 
Js e (l) 

This integrand is holomorphic near u = 1 for to < n , and the integral therefore vanishes in that case. For to > n we 
use the invariance of the integration form under u — > —u to write I m ,n as an integral over a sum of two circles: 

W = Sri I ( U 2 = l )m -n + l • 1 = S ^ + • 

The integrand in this case is holomorphic near u = . In the punctured plane C \ {1, —1} the chain S e (l) + S e (— 1) 
is homologous to the circle at infinity, where the integrand vanishes. Therefore the integral again is zero. This proves 
that I m ,n = for to ^ n, which in turn completes our proof that the contour integrals (I5.20|) and l|5.21|l are the same 
as the functions P n (uj), Q n {oj) defined from 15.71 15. 8|) by (|5.12l I5?T3|) . As a final check, note that 

In,n = (TTir 1 / / f 1 = 1 , 

/s e (i) u ( u2 - !) 

which is what it ought to be in view of l)5.14[l . 

Now we harvest a major benefit from the contour integral representations l|5.20|l and (|5.21|) : using these, we can 
carry out the sum in the definition 15.16f) of the kernel Kn as a geometric sum. The result is a double contour integral: 

Kn(uji,uj2) = <j> du <j> dv Fjy(u,v ;uji 7 u}2) , (5.25) 

S e (l) SAO) 

This exact expression represents the complete solution of our problem. We will now use it to determine the large- N 
asymptotics in the bulk and at the hard edge u> = . 



D. Asymptotics in the bulk 

The kernel on the diagonal u>\ = uj-2 is the same as the 1-level function, Ri(u) = Kn(uj,uj); sec l|5.19ll . We already 
know from Sect. II Vl the asymptotics of Ri(to) = p(co) in the bulk: introducing the scaling variable x = lo/N (formerly 
x = ut/N), this is 

lim K N (Nx,Nx) = Poc(x) , 

with Poo(x) given by l|4.6|l . In the present subsection we are going to demonstrate that the scaling limit of the kernel 
Kn(lui,lu2) off the diagonal leads to sine-kernel universality for all level correlation functions: 

lim R n (Nx + w 1 ,...,Nx + u n ) = Det v , - \ — ^- , (5.27) 

N^oc \ 7r(o>i — LUj) I 

\ / — l,...,n 

as is expected for systems in the universality class of the Gaussian Unitary Ensemble. As a corollary, we will obtain 
an independent confirmation of the result l|4.6() . 

Looking at the integral representation l|5.25|l one might think that the large- N limit could be taken by applying 
the saddle-point method to that integral. However, as we shall sec, the dominant saddle points lie on the line u = v 
where the integrand has a singularity of type 0/0 which, albeit removable, complicates the saddle-point evaluation. 
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Therefore, rather than calculating K^(uj\, 0J2) directly, we look at the product (uji — L02) Kn(uii, uj2)- Using the 
relation (u 2 — ^i) e" 21 '^ 1 '" = (d v + d u ) e^ 2 "-^" an d partially integrating, we rewrite l|5.25f) as 

(oji - uj 2 ) K n {uji,uj2) — j> du j> dv F N (u,v ;wi,w 2 ) , (5.28) 

S e (l) S e (0) 

1 ,, .. ( o a \ u i v- i+i ( n-v-° 



F »^; W) = ^- M -\tu + f v )¥^{{r^) -*)> (5 ' 29) 

which constitutes the starting point for the following analysis. 

In preparation for taking the limit TV — > 00 , we set ui\ = Nx + uj and uj 2 = Nx + Gj . The deciding factor in the 
integrand of [w\ — w 2 ) Kn(u)\ , tt> 2 ) m the large- TV limit will then be 

exp (-Nx{u-v) + TVlog(l - v~ 2 ) - TV log (1 - u~ 2 )) , 

leading to the saddle-point equation 

—2 

ip{u) = x = <p(v) , ip{w) = -9^1og(l - w~ 2 ) = — — . (5.30) 

w(w z — 1) 

Notice that ip is related to our function i|4.17[) by f(w — 1) = ip(w). A comprehensive study of the equation f(w) = z 
and its solutions for w was made in Sect. llVI From there we know that the saddle-point equation <p(w) = x has three 
solutions in general, and for < x < b = 3\/3 these are 

w a (x) = -x-^ 3 e- 2nia ^{l + ^/l-x 2 /b 2 ) 1/3 - x- 1/3 e 27ria/3 (l - ^/l-x 2 /b 2 ) 1/3 (a = 1, -1, 0) . (5.31) 

In the range of interest (0 < x < b) the first two solutions, w±i(x), are complex conjugates of each other while the 
third one, wo(x), is negative. Expanding the logarithm of (1 — v~ 2 ) N /(l — u~ 2 ) N to second order around a pair of 
saddle points w a , wp we encounter the Gaussian 

f>7/J 2 — 9 

exp (±N<p'(w a (x))(6u) 2 - \N^{w B {x))(8v) 2 ) , <p'(w) 



w 2 (w 2 — l) 2 

For the negative saddle point one has <p'(wo(x)) > , so its path of steepest descent would be perpendicular to the 
real axis in the case of u and along the real axis in the case of v. The latter is inconsistent with the original integration 
contour for v being S c (0). In the former case, wq(x) < —1 is inaccessible because of the singularity of (1 — u~ 2 )~ N 
intervening at it = — 1. Thus this saddle point is irrelevant for present purposes and may be discarded. 

We now make another preparation of the saddle-point evaluation of the integral, by investigating the behavior of 
the integrand near the two remaining saddle points. We set 

u = w a (x) + N~ 1/2 8u , v = wp(x) + N~ 1/2 8v (a,/3 = ±l), 

and first look at the diagonal case, a = (3 . Using the identity 

1 (8 d\(\-v- 2 \ N n ^ T u 2 +uv + v 2 -l (l-v- 2 \ N , oo , 

hr + T" 1 =2 = 2N t 2 1W 2 TT i =a ' 5 - 32 

u — v \ou ov J \1 — u z J u(u z — 1) v(v z — 1) \ 1 — u 1 J 

we find the scaling limit of the integrand Fn to be 

lim N^FNiwaix) + N- 1/2 Su,w Q (x) + N~ 1/2 5v ; Nx + lj , Nx + u) 

TV— »OC 

= (2ny 2 e- w ^ x) ^-^ip\w a (x))c^' {Wa{x)]{Su2 - 5v2) . 
The same limit in the off-diagonal case (a ^ (3) vanishes. Indeed, 

2^ .2 -1 w a (wl-l)-w p (w}-l) -2 _ x _! 
w a + w a WB + Wg - 1 = = (x -x ) = laf^p), 

W a - Wp W a - WB 

and therefore the factor in the numerator on the right-hand side of 1)5. 32|) gives zero. 
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FIG. 2: Sketch of the saddle points for the case of x = 1 . By deforming the original contours, which are small circles around 
the singular points u = 1 and v — 0, one arranges for the contours of integration to pass through the saddle points w±i(x) in 
the direction of steepest descent. Away from the saddle points the deformed contours are drawn as paths of constant phase, 
which interpolate between different zeroes of the integrand: they run between and +oo for u , and between 1 and -co for v . 

We now deform the contours of integration as indicated in Fig. [5J The deformed contours pass through the saddle 
points w±i but miss the saddle point wq ■ At w±i the paths of steepest descent for u and v cross at right angles, 
valleys in one case being mountains in the other case and vice versa. 

Next we do the Gaussian integrals. Given the counterclockwise orientations of the original contours S e (l) resp. 
S e (0), and taking into account the directions of the paths of steepest descent, we get 

The product of these two integrals is 2-k\/ ip 1 (w+i). The same calculation for the other saddle w-\ gives — 2ivi/ip'(w-i). 
Thus, putting the factors together and summing over the contributions from diagonal pairs of saddle points (a = (3) 
we obtain 

(u-w) lim K n (Nx + cj,Nx + u>) = — ( fl -«-iWM_ fl -«nWM)) . 

Since [Heiu+i = D\tw~i and 3mw + i = — 3mw_i this means that, with 3mw+i(x) =: Trp 00 (x), we have 

ns i \r -\ sin (■kp„ Ci {x)(u} — £))) 
lim K N (Nx + u,Nx + Q) = c -*"»±i(*)(»-") x v H ™ K )K . '± . 

N—*oo Tt{u> — Uj) 

The exponential factor e _5Re w ±i( x )( UJ - UJ ) drops out when forming the determinant on the right-hand side of l|5.19|l . 
Thus we arrive at the universal sine-kernel (or GUE) correlation functions (|5.27(l . 

Setting lj — uj notice the special result limjy^oo Kn(Nx,Nx) = Poo(x) . Since the kernel on the diagonal is none 
other than the 1-level function, Kn(lu,oj) = p(u>), this gives another determination of the large-iV level density poo . 
From l|5.31(l one sees that poo{x) = w~ 1 3mw+i(x) agrees with our earlier result l|4.ri|l . 

E. Asymptotics near uj — 

At the lower edge (uj = 0) of the spectrum, a new type of behavior is expected to emerge. This behavior, as we 
shall see presently, occurs on a scale uj ~ N" 1 ^ 2 . 



^\<p'{w+i)\- 



l / 2 e l(- 1ri - iar S l P'( u '+i)) 



^\<p'(w +1 )\- 



1 / 2 e ^(27ri-iarg ip' (w +1 )) 
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y y 
FIG. 3: The graph of the scaling function k(y,y) for the case of I — (left) and I = 1 (right). 



To exhibit the scaling limit near co = , it is best to send the integration variables u, v to their reciprocals, u — > u^ 1 
and v — > v -1 . Then c?u — > ~u~ 2 du, dv — > —v~ 2 dv, and the integration contour for u has its radius inverted and 
orientation reversed, 5 e (0) — ♦ — Si/ e (0). However, since the integrand is holomorphic in uonC\ {0} we may return 
to the original radius e (or any other radius, for that matter). In the case of u we take the radius e of 5 e (l) to be 
very small. Then inversion u — > w _1 sends S^l) to itself (or, in any case, to the same homology class on C \ {1}), 
with no change of orientation. Altogether, then, carrying out the transformation (u,v) — ► (u ,v ) the integral 
representation 1|5.25[) continues to hold true if we make the replacement 



F n (u,v;uj 1 ,u 2 ) 



2 v 2 F N {u l ,v 1 ;ui,u 2 ) 



1 p -Wl/«+UJ2/u U l V l 1 



2tt 2 



1 - W 



N 



Next, as another preparation for taking the limit TV — > oo, we deform the u-contour ^(l) to some axis parallel 
to the imaginary axis. The deformed contour crosses the real axis between u = and u = 1 and is directed from 
u = +ioo to u = — ioo. We also reverse the direction of integration for u and change the overall sign of the integral. 

Then we set Uj = N~ x l 2 yj and rescale u — > N~ 1 / 2 u and v — > N~ 1 / 2 v accordingly. Again, in view of the analytic 
properties of the integrand we can keep the integration contours fixed while rescaling. Because the u-integral converges 
at infinity we have a good limit 



lim (1 - u /N) 



-N 



exp(u 2 ) 



In total, we thus obtain the following scaling limit for our kernel Kn ■ 

1 f , f dv 



k(yi,y 2 ):= lim N- 1 / 2 K N {N- 1 / 2 y 1 ,N- 1 ' 2 y 2 ) 

iV— i-oc 



2, 2 J dU f~ C 

iR+e Ui 



-yi/u+y 2 /v 



(v/u) 1 



(5.33) 



where Ui = £"1(0) means the unitary numbers, and iR + e is the imaginary axis translated by e > into the right 
half of the complex plane. Plots of the scaling function k(y,y) for I = 0, 1 arc shown in Fig. 01 Using the method of 
saddle-point evaluation as in Sect. IV 01 one can show that this function behaves as k(y,y) ~ y~ 1 ^ 3 for large y. 
Taking the same scaling limit for the functions P/v(w) and Qn(w) in 15.2l)f) and (|5.21|) one gets 



p{y) = lim N~( l ~ 1 ^ 2 P N (N~ 1 / 2 y) 



1 

?ri 



,u 2 -y/u u -l du 



(5.34) 



q(y) = lim N l / 2 Q N (N~y 2 y) 



Both functions have convergent series expansions: 

(-yy 



p 



(y) = E 



n\T((l + n + l)/2) ' 



2?ri 



e -v 2 +y/v v l-l dv 



Ui 



i-y 2 ) 



2\n 



n=0 



n\ (2n + l)\ 



(5.35) 



(5.36) 
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The expansion for q(y) can be obtained either directly from (|5.35(l . or by taking the limit N — > oo in (|5.22() . In 
the case of p(y), the earlier formula (|5.23[) is not suitable; rather, in order to verify (|5 . 3(jp> for p(y) one expands the 
integrand of (|5.34(1 in powers of y , makes use of 9\e u = t > to write 



does the Gaussian it-integral by completing the square, and uses the duplication formula for the Gamma function. 

Acknowledgement. We thank our colleague T. Kricchcrbaucr for helpful discussions. This work has been supported 
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Note added. After submission of this manuscript, P. Forrester pointed out to us that the joint eigenvalue distribution 
derived and analyzed here falls in a broad class of models solved by Borodin ^3 ■ Borodin's expression for the kernel 
K]y(x,y) is equivalent to ours by old work of Konhauser |l3l |. The mathematical results of Konhauser were first 
introduced into random matrix physics by Muttalib |l4j , who suggested to use them for an approximate treatment of 
the statistics of transmission eigenvalues of disordered conductors. 
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